######################################################################
#Replication file for "The Impact of China's AIIB on the World Bank"
#Authors: Jing Qian, James Vreeland, Jianzhi Zhao
#Codes to replicate results in Appendix A, column (5), and Figure A.1
######################################################################
rm(list = ls())

#To install the bpCausal package (ver 0.0.1)
# install.packages('devtools', repos = 'http://cran.us.r-project.org') # if not already installed
# devtools::install_github('liulch/bpCausal')

#gsynth package can be installed from CRAN (ver 1.1.8)
# install.packages("gsynth")


#Load bpCausal package
library(bpCausal)

#Load gsynth package
library(gsynth)

#Other packages used
library(tidyverse)

#Load custom functions
#Note: Given the number of tests performed in the appendix, 
#several custom functions are created to implement estimation methods in a more efficient way.
#These functions merely create a wrap-up of existing functions, like bpCausal and gsynth,
#in order to avoid typing parameters repeatedly.

source("code/custom_functions.R")

#Setwd
setwd("C:/Users/qianj/Dropbox (Princeton)/AIIB_WB_Replication")


###################################
#Table A.1, column (5) 
#3 minutes
###################################
#-------------------------
#(1) Load data
#-------------------------
load("data/data_main.RData")

#-------------------------
#(2) Parameters
#-------------------------
iv.use = "aiib_founder_2016"
dv.list = c("hard_project_count_all",
            "soft_project_count_all")
covs = c("gdppc_log_lag",
         "population_log_lag",
         "election_either_lag",
         "fdi_gdp_lag",
         "debt_gni_lag",
         "oda_gni_lag",
         "polity2_lag",
         "unsc_lag",
         "IdealPointDistance_lag")

#gsynth requires complete dataset
df.use = df.main[, c(dv.list, iv.use, covs, c("iso2c", "year"))] %>%
  .[complete.cases(.), ]

#-------------------------
#(3) Estimation
#-------------------------
# result.gsynth = gsynth.group(dv.list = dv.list,
#                              iv = iv.use,
#                              covs = covs,
#                              df = df.use)


result.gsynth = list()
result.gsynth[[1]] = gsynth(Y = dv.list[1],
                            D = iv.use,
                            X = covs,
                            data = df.use,
                            index = c("iso2c", "year"),
                            r = c(0, 10),
                            force = "two-way",
                            se = TRUE,
                            nboots = 2000,
                            CV = TRUE,
                            inference = "parametric",
                            seed = 1234,
                            min.T0 = 12)
result.gsynth[[2]] = gsynth(Y = dv.list[2],
                            D = iv.use,
                            X = covs,
                            data = df.use,
                            index = c("iso2c", "year"),
                            r = c(0, 10),
                            force = "two-way",
                            se = TRUE,
                            nboots = 2000,
                            CV = TRUE,
                            inference = "parametric",
                            seed = 1234,
                            min.T0 = 12)
names(result.gsynth) = dv.list

#-------------------------
#(4) Estimated coefficients and intervals
#-------------------------
att.tab.a1.5 = result.gsynth$hard_project_count_all$est.avg #Estimated ATT as reported in Table A.5, column (5)
att.tab.a1.5 %>%
  round(digits = 3) %>%
  format(nsmall = 3) #Rounding, as reported in Table A.1, column (5)

coef.tab.a1.5 = result.gsynth$hard_project_count_all$est.beta #Estimated covariate coefficients as reported in Table A.5, column (5)
coef.tab.a1.5 %>%
  round(digits = 3) %>%
  format(nsmall = 3) #Rounding, as reported in Table A.1, column (5)

#-------------------------
#(5) Number of observations
#-------------------------
#Three units removed due to too few pre-treatment periods: AF, DJ, TL
df.use.complete = df.use %>%
  subset(!iso2c %in% c("AF", "DJ", "TL"))
nrow(df.use.complete)

#-------------
#(6) Treated and Control Units
#-------------
length(result.gsynth$hard_project_count_all$id.tr) #Treated units
length(result.gsynth$hard_project_count_all$id.co) #Control units

#-------------
#(7) Unobserved factors
#-------------
result.gsynth$hard_project_count_all$r.cv


###################################
#Figure A.1, Estimate using bpCausal
###################################

#-------------
#(1) Initialization
#-------------
load("data/data_main.RData") #Load data
dv.list = c("hard_project_count_all",
            "soft_project_count_all")
D.use = "aiib_founder_2016"
X.use = Z.use = A.use = c("gdppc_log_lag",
                          "population_log_lag",
                          "election_either_lag",
                          "fdi_gdp_lag",
                          "debt_gni_lag",
                          "oda_gni_lag",
                          "polity2_lag",
                          "unsc_lag",
                          "IdealPointDistance_lag")

#-------------
#(2) Estimation
#5.4 minutes
#-------------
result.bpcausal = bpcausal.group(dv.list = dv.list,
                                 df.use = df.use,
                                 D.use = D.use,
                                 X.use = X.use,
                                 Z.use = Z.use,
                                 A.use = A.use)

###################################
#Figure A.1
###################################
#--------------------
#(1) Get results
#--------------------
result.hard.use = effSummary(result.bpcausal$hard_project_count_all)
result.soft.use = effSummary(result.bpcausal$soft_project_count_all)

gsynth.hard.use = result.gsynth$hard_project_count_all
gsynth.soft.use = result.gsynth$soft_project_count_all

#-------------
#(2) Parameters
#-------------
#Titles
xtitle = ""
ytext = ""

#Limits
ylim = c(-3, 3)
newx = 1992:2019
t.time = 2016

#Cex
cex = 2.2


#-------------
#(4) Plot
#-------------
#Setup
png(filename = "figure/Figure_A1.png",
    height = 1600,
    width = 2400)

par(mfrow = c(2, 2),
    mar = c(5, 5, 5, 1))

################
#bpCausal, infrastructure
################
#----------
#Plot empty graph
#----------
plot(1,
     type = "n",
     xlab = "",
     ylab = "",
     axes = F,
     xlim = range(newx),
     ylim = ylim)
box()

title(main = "DM-LFM: World Bank Infrastructure Projects",
      line = 1,
      cex.main = cex)

#-----------------------
#Add axis
#-----------------------
#X-axis labels
axis(side = 1, 
     at = seq(min(newx),
              max(newx),
              by = 2),
     cex.axis = cex)
#X-axis title
mtext(text = xtitle, 
      side = 1, 
      line = 2.5,
      cex = cex)
#Y-axis
axis(side = 2,
     at = seq(min(ylim),
              max(ylim),
              by = 1),
     cex.axis = cex)
#Y-axis title
mtext(text = ytext,
      side = 2,
      line = 2.5, 
      cex = cex)

#-----------------------
#Lines and shades
#-----------------------
#Ablines
abline(v = t.time,
       col = "gray",
       lty = 2,
       lwd = cex)
abline(h = 0,
       col = "gray20",
       lty = 2,
       lwd = cex)

#ATTs
lines(x = newx,
      y = result.hard.use$est.eff$estimated_ATT,
      col = 1,
      lty = 1,
      lwd = cex)
#CI
polygon(x = c(rev(newx), 
              newx),
        y = c(rev(result.hard.use$est.eff$estimated_ATT_ci_l), 
              result.hard.use$est.eff$estimated_ATT_ci_u),
        col = "#55555530",
        border = NA)
#Legend
legend("topleft",
       legend = c("Estimated ATT",
                  "95% Credible Intervals"),
       cex = cex,
       seg.len = 2,
       col = c(1, "#55555530"),
       lty = c(1, 5),
       lwd = c(2, 20),
       bty = "n")


################
#gsynth, infrastructure
################
#----------
#Plot empty graph
#----------
plot(1,
     type = "n",
     xlab = "",
     ylab = "",
     axes = F,
     xlim = range(newx),
     ylim = ylim)
box()

title(main = "GSC: World Bank Infrastructure Projects",
      line = 1,
      cex.main = cex)

#-----------------------
#Add axis
#-----------------------
#X-axis labels
axis(side = 1, 
     at = seq(min(newx),
              max(newx),
              by = 2),
     cex.axis = cex)
#X-axis title
mtext(text = xtitle, 
      side = 1, 
      line = 2.5,
      cex = cex)
#Y-axis
axis(side = 2,
     at = seq(min(ylim),
              max(ylim),
              by = 1),
     cex.axis = cex)
#Y-axis title
mtext(text = ytext,
      side = 2,
      line = 2.5, 
      cex = cex)

#-----------------------
#Lines and shades
#-----------------------
#Ablines
abline(v = t.time,
       col = "gray",
       lty = 2,
       lwd = cex)
abline(h = 0,
       col = "gray20",
       lty = 2,
       lwd = cex)

#ATTs
lines(x = newx,
      y = gsynth.hard.use$est.att[, 1],
      col = 1,
      lty = 1,
      lwd = cex)
#CI
polygon(x = c(rev(newx), 
              newx),
        y = c(rev(gsynth.hard.use$est.att[, 3]), 
              gsynth.hard.use$est.att[, 4]),
        col = "#55555530",
        border = NA)
#Legend
legend("topleft",
       legend = c("Estimated ATT",
                  "95% Confidence Intervals"),
       cex = cex,
       seg.len = 2,
       col = c(1, "#55555530"),
       lty = c(1, 5),
       lwd = c(2, 20),
       bty = "n")

################
#bpCausal, non-infrastructure
################
#----------
#Plot empty graph
#----------
plot(1,
     type = "n",
     xlab = "",
     ylab = "",
     axes = F,
     xlim = range(newx),
     ylim = ylim)
box()

title(main = "DM-LFM: World Bank Non-Infrastructure Projects",
      line = 1,
      cex.main = cex)

#-----------------------
#Add axis
#-----------------------
#X-axis labels
axis(side = 1, 
     at = seq(min(newx),
              max(newx),
              by = 2),
     cex.axis = cex)
#X-axis title
mtext(text = xtitle, 
      side = 1, 
      line = 2.5,
      cex = cex)
#Y-axis
axis(side = 2,
     at = seq(min(ylim),
              max(ylim),
              by = 1),
     cex.axis = cex)
#Y-axis title
mtext(text = ytext,
      side = 2,
      line = 2.5, 
      cex = cex)

#-----------------------
#Lines and shades
#-----------------------
#Ablines
abline(v = t.time,
       col = "gray",
       lty = 2,
       lwd = cex)
abline(h = 0,
       col = "gray20",
       lty = 2,
       lwd = cex)

#ATTs
lines(x = newx,
      y = result.soft.use$est.eff$estimated_ATT,
      col = 1,
      lty = 1,
      lwd = cex)
#CI
polygon(x = c(rev(newx), 
              newx),
        y = c(rev(result.soft.use$est.eff$estimated_ATT_ci_l), 
              result.soft.use$est.eff$estimated_ATT_ci_u),
        col = "#55555530",
        border = NA)
#Legend
legend("topleft",
       legend = c("Estimated ATT",
                  "95% Credible Intervals"),
       cex = cex,
       seg.len = 2,
       col = c(1, "#55555530"),
       lty = c(1, 5),
       lwd = c(2, 20),
       bty = "n")

################
#gsynth, non-infrastructure
################
#----------
#Plot empty graph
#----------
plot(1,
     type = "n",
     xlab = "",
     ylab = "",
     axes = F,
     xlim = range(newx),
     ylim = ylim)
box()

title(main = "GSC: World Bank Non-Infrastructure Projects",
      line = 1,
      cex.main = cex)

#-----------------------
#Add axis
#-----------------------
#X-axis labels
axis(side = 1, 
     at = seq(min(newx),
              max(newx),
              by = 2),
     cex.axis = cex)
#X-axis title
mtext(text = xtitle, 
      side = 1, 
      line = 2.5,
      cex = cex)
#Y-axis
axis(side = 2,
     at = seq(min(ylim),
              max(ylim),
              by = 1),
     cex.axis = cex)
#Y-axis title
mtext(text = ytext,
      side = 2,
      line = 2.5, 
      cex = cex)

#-----------------------
#Lines and shades
#-----------------------
#Ablines
abline(v = t.time,
       col = "gray",
       lty = 2,
       lwd = cex)
abline(h = 0,
       col = "gray20",
       lty = 2,
       lwd = cex)

#ATTs
lines(x = newx,
      y = gsynth.soft.use$est.att[, 1],
      col = 1,
      lty = 1,
      lwd = cex)
#CI
polygon(x = c(rev(newx), 
              newx),
        y = c(rev(gsynth.soft.use$est.att[, 3]), 
              gsynth.soft.use$est.att[, 4]),
        col = "#55555530",
        border = NA)
#Legend
legend("topleft",
       legend = c("Estimated ATT",
                  "95% Confidence Intervals"),
       cex = cex,
       seg.len = 2,
       col = c(1, "#55555530"),
       lty = c(1, 5),
       lwd = c(2, 20),
       bty = "n")


#----------
#close
dev.off()
